/***********************************************************************/
/* Open Visualization Data Explorer                                    */
/* (C) Copyright IBM Corp. 1989,1999                                   */
/* ALL RIGHTS RESERVED                                                 */
/* This code licensed under the                                        */
/*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
/***********************************************************************/

#include <dxconfig.h>

#include <string.h>
#include <stdarg.h>
#include "arrayClass.h"

ProductArray _NewProductArrayV( int n, Array *terms,
                                struct productarray_class *class )
{
  Type type;
  Category category;
  ProductArray a;
  int i, items, rank, *shape;

  /* set type according to first term */
  if ( n > 0 )
  {
    type = terms[0]->type;
    category = terms[0]->category;
    rank = terms[0]->rank;
    shape = terms[0]->shape;
  }
  else
  {
    type = TYPE_FLOAT;
    category = CATEGORY_REAL;
    rank = 0;
    shape = NULL;
  }

  /* create array */
  a = (ProductArray)_dxf_NewArrayV( type, category, rank, shape,
                                    (struct array_class *)class );
  if ( !a )
    return NULL;

  /* typecheck and set terms */
  for ( i = 0, items = 1; i < n; i++ )
  {
    if ( !DXTypeCheckV( terms[i], type, category, rank, shape ) )
      DXErrorReturn( ERROR_BAD_TYPE, "Inconsistent product array term types" );
    items *= terms[i]->items;
  }

  /* copy info */
  a->nterms = n;
  if ( n * sizeof( *terms ) <= sizeof( a->array.ldata ) )
  {
    a->terms = (Array *)a->array.ldata;
  }
  else
  {
    a->terms = (Array *)DXAllocate( n * sizeof( *terms ) );
    if ( !a->terms )
      return NULL;
  }
  for ( i = 0; i < n; i++ )
    a->terms[i] = (Array)DXReference( (Object)terms[i] );
  a->array.items = items;

  return a;
}

ProductArray DXNewProductArrayV( int n, Array *terms )
{
  return _NewProductArrayV( n, terms, &_dxdproductarray_class );
}

ProductArray DXNewProductArray( int n, ... )
{
  Array terms[100];
  int i;
  va_list arg;

  ASSERT( n < 100 );

  /* collect args */
  va_start( arg, n );
  for ( i = 0; i < n; i++ )
    terms[i] = va_arg( arg, Array );
  va_end( arg );

  /* call V version */
  return DXNewProductArrayV( n, terms );
}

ProductArray DXGetProductArrayInfo( ProductArray a, int *n, Array *terms )
{
  int i;

  CHECKARRAY( a, CLASS_PRODUCTARRAY );

  if ( n )
    *n = a->nterms;

  if ( terms )
    for ( i = 0; i < a->nterms; i++ )
      terms[i] = a->terms[i];

  return a;
}

/*
 * XXX - the DXsyncmem()s below are because somebody else may have created
 * the data and the allocation for it might be in a new segment.  This
 * is little kludgy; we really need a better solution to this problem.
 */

Pointer _dxfProductArray_GetArrayData( ProductArray a )
{
  int i, j, k, l, nterms, ndim, nl, nr, rank;
  Array *terms;
  float *result = NULL, *res, *right, *rt, *left, *lt;
  int freelocal = 0;

  /* DXlock array */
  EXPAND_LOCK( a );

  nterms = a->nterms;
  terms = a->terms;
  rank = a->array.rank;
  for ( i = 0, ndim = 1; i < rank; i++ )
    ndim *= a->array.shape[0];

  /* first left operand is first term */
  left = (float *)DXGetArrayData( terms[0] );
  if ( !left )
    goto error;
  nl = terms[0]->items;
  /*ASSERT(terms[0]->shape[0]==ndim);*/

  /* copy first term data if nterms is 1 */
  /* XXX - is there any way to avoid this? */
  if ( nterms == 1 )
  {

    int size = nl * ndim * sizeof( float );
    result = (float *)DXAllocate( size );
    if ( !result )
      goto error;
    memcpy( result, left, size );

    /* accumulate product left to right */
  }
  else
    for ( i = 1; i < nterms; i++ )
    {

      /* right operand is next term */
      right = (float *)DXGetArrayDataLocal( terms[i] );
      if ( !right )
      {
        freelocal = 0;
        right = (float *)DXGetArrayData( terms[i] );
        if ( !right )
          goto error;
      }
      else
        freelocal = 1;

      nr = terms[i]->items;
      /*ASSERT(terms[i]->shape[0]==ndim);*/

      /* allocate result */
      res = result = (float *)DXAllocate( nl * nr * ndim * sizeof( float ) );
      if ( !result )
        goto error;

      /* left operand plus right operand */
      switch ( ndim )
      {
        case 1:
          for ( j = 0, lt = left; j < nl; j++, lt += ndim )
          {/* ea left point */
            float t0 = lt[0];
            for ( k = 0, rt = right; k < nr;
                  k++, rt += ndim ) /* ea right point */
              *res++ = t0 + rt[0];  /* sum */
          }
          break;

        case 2:
          /* left operand plus right operand */
          for ( j = 0, lt = left; j < nl; j++, lt += ndim )
          {/* ea left point */
            float t0 = lt[0];
            float t1 = lt[1];
            for ( k = 0, rt = right; k < nr; k++, rt += ndim )
            {                      /* ea right point */
              *res++ = t0 + rt[0]; /* sum */
              *res++ = t1 + rt[1]; /* sum */
            }
          }
          break;

        case 3:
          /* left operand plus right operand */
          for ( j = 0, lt = left; j < nl; j++, lt += ndim )
          {/* ea left point */
            float t0 = lt[0];
            float t1 = lt[1];
            float t2 = lt[2];
            for ( k = 0, rt = right; k < nr; k++, rt += ndim )
            {                      /* ea right point */
              *res++ = t0 + rt[0]; /* sum */
              *res++ = t1 + rt[1]; /* sum */
              *res++ = t2 + rt[2]; /* sum */
            }
          }
          break;

        default:
          /* left operand plus right operand */
          for ( j = 0, lt = left; j < nl; j++, lt += ndim )
          {/* ea left point */
            for ( k = 0, rt = right; k < nr;
                  k++, rt += ndim )        /* ea right point */
              for ( l = 0; l < ndim; l++ ) /* ea dimension */
                *res++ = lt[l] + rt[l];    /* sum */
          }
          break;
      }

      /* next left operand is result */
      if ( i > 1 )
        DXFree( (Pointer)left );
      left = result;
      nl *= nr;
      if ( freelocal )
        DXFreeArrayDataLocal( terms[i], (Pointer)right );
    }

  ASSERT( nl == a->array.items );

  /* DXunlock and return */
  EXPAND_RETURN( a, result );

error:
  if ( freelocal )
    DXFreeArrayDataLocal( terms[i], (Pointer)right ); /* and fall thru */
  EXPAND_ERROR( a );
}

Error _dxfProductArray_Delete( ProductArray a )
{
  int i;
  for ( i = 0; i < a->nterms; i++ )
    DXDelete( (Object)a->terms[i] );
  if ( (Pointer)a->terms != (Pointer)a->array.ldata )
    DXFree( (Pointer)a->terms );
  return _dxfArray_Delete( (Array)a );
}
